#include <NTL/ZZ_pX.h>

/*
Compile with: 
$ ~/Desktop/Summer2010/Poly_Factor$ g++  -I$HOME/sw/include test.cpp -o foo  -L$HOME/sw/lib -lntl  -lm
*/

NTL_CLIENT

void DDF(ZZ_pX,long);
long EDF(ZZ_pX,long,long);

int main()
{
	long p;
	cout << "Enter Radix: ";
	cin >> p;
	ZZ radix;
	radix = p;
	ZZ_p::init(radix);

	ZZ_pX f;
	cout << "Enter Polynomial module " << p << " :";
	cin >> f;

	cout << "Polynimial: " << f << endl;	
#if 0
	cout << "Calling DDF ... " << endl;
	DDF(f,p);
	cout << "Returning from DDF ... " << endl;
#endif
	cout << "Calling EDF ..." << endl;
	EDF(f,p,1);
	cout << "Returning from EDF ..." << endl;

	//vec_pair_ZZ_pX_long factors;

	//CanZass(factors, f);  // calls "Cantor/Zassenhaus" algorithm

	//cout << factors << "\n";
    
}

/**
 * Distinct Degree Factorization
 */
void DDF(ZZ_pX f,long p)
{
	ZZ_pX x; //x = X
	SetX(x);

	ZZ_pX h_0 = power(x,p);	// h_0 = x^p 

	ZZ_pX X_eq_1;
	set(X_eq_1);	// X_eq_1 = 1

	long i;
	long degree_of_f = deg(f);

	for(i=1; i < degree_of_f; i=i+1)
	{
		ZZ_pX h = power(h_0,i);
		//cout << "h_" << i << " : " << h << endl;
		ZZ_pX g = GCD(f,h - x);	// g = gcd(f,(x^p)^i - x)
		if( g != X_eq_1 )
		{
			cout << "g_" << i << " :" << g << endl;
			f = f/g;
		}
	}
}
/**
 * Equal Degree Factorization
 */
long EDF(ZZ_pX a,long p,long n)
{
cout << "EDF called on: " << a << endl;

	long num_of_factors = 0;

	if( deg(a) <= n )	
	{
		num_of_factors++;
		cout << "***Factor: " << a << endl;
cout << "# of factors found: " << num_of_factors << endl;
		return(num_of_factors);
	}
	long	m = deg(a)/n;	// m  = number of factors of degree n in a(x)
cout << endl << "m: " << m << endl << endl;

	// construct 1
	ZZ_pX X_eq_1;
	set(X_eq_1);	// X_eq_1 = 1

	while(num_of_factors < m)
	{
		long rand_deg = rand() % (2*n);	// rand_deg = the degree of random polynomial v(x) 		
		ZZ_pX v = random_ZZ_pX(rand_deg+1);	// v(x) = random polynomial of degree at most (2*i-1)
cout << "Random Poly of degree *" << rand_deg << "*: " << v << endl;

		// if p is odd
		v = power(v, (p^n-1)/2 )-1; // v(x) = v(x)^((p^n-1)/2) - 1;
cout << "v(x)^((p^n-1)/2): " << v << endl; 
		ZZ_pX g = GCD(a,v);
cout << "gcd with *" << a << "*: " << g << endl;

		if( (g != X_eq_1) && (g != a) )
		{
cout << "CHECK POINT 1" << endl;
int CP1; cin >> CP1;
			long temp1 = EDF(g,p,n); 
			long temp2 = EDF(a/g,p,n);			
			num_of_factors = temp1 + temp2; 
cout << "## of factors found:" << num_of_factors << endl;
cin >> CP1;
				
		}
	}

	return num_of_factors;
}

